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We give a comprehensive introduction into a diagram¬ 
matic method that allows for the evaluation of Gutzwiller 
wave functions in finite spatial dimensions. We discuss 
in detail some numerical schemes that turned out to be 
useful in the real-space evaluation of the diagrams. The 
method is applied to the problem of d-wave 


superconductivity in a two-dimensional single-band 
Hubbard model. Here, we discuss in particular the role 
of long-range contributions in our diagrammatic expan¬ 
sion. We further reconsider our previous analysis on the 
kinetic energy gain in the superconducting state. 
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1 Introduction For the study of the ground-state 
properties of quantum systems, variational wave func¬ 
tions can be a powerful tool. The most common variational 
approach in the theory of correlated electron systems is 
the Hartree-Fock approximation (HFA) which is based on 
(variational) single-particle product wave functions |<fo)- If 
applied to systems with attractive, e.g., phonon-mediated 
interactions, the HF theory leads to the celebrated BCS 
theory on superconductivity ID. 

The crucial first step in any variational approach is 
the calculation of the energy expectation value for a given 
class of variational wave functions. For HF wave functions 
this can always be achieved by means of Wick’s theorem, 
which explains the popularity of this approach. Many phe¬ 
nomena in correlated electron systems, however, cannot be 
described properly by the HFA, especially, in (effectively) 
one- or two-dimensional systems. This holds, in particular, 
for unconventional superconductivity which has been ob¬ 
served in a number of materials, such as Cuprates, Ruthen- 
ates and iron-based Pnictides. 

A way to improve the HFA is based on ‘Jastrow wave 
functions’ which have the form l2ll3 

I'f'j) = Al'I'o) . (1) 


Here, Pj is an operator which has been chosen in various 
ways in the literature IHIIsIIMTIISI IWTOII and is meant to ac¬ 
count for correlation effects which are not captured by the 
Hartree-Fock (single-particle) wave function jifb). One of 
the simplest examples for such a Jastrow wave function is 
the Gutzwiller wave function ||4ll5ll^ which will be used in 
this work. 

Evaluating expectation values for Jastrow (or Gutz¬ 
willer) wave functions is a difficult many-particle prob¬ 
lem which, in general, can only by tackled by numerical 
techniques, such as the ‘variational Monte-Carlo method’ 
(VMC) ifTnfTHfra . We have recently developed a dia¬ 
grammatic scheme for the evaluation of expectation val¬ 
ues for Gutzwiller wave functions. Unlike the VMC, our 
method adresses the infinite systems and, hence, it does 
not suffer from the typical finite-size errors of VMC. As 
shown in Ref. na , our approach allows us to study e.g., 
the stability of nematic (‘Pomeranchuk’) phases in two- 
dimensional Hubbard models. First results on the stabil¬ 
ity of superconducting ground states in these models have 
been presented in Ref. ifTSll . 

In this work we will give a comprehensive introduction 
into the technical details of our approach for the study of 
superconducting ground states and present numerical re- 
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suits which complement those published in previous work. 
Refs. i[ii[i5i[i6i[ni[ia . We will introduce, in particular, a 
new way to evaluate diagrams which contain long-range 
correlations. 

Our presentation is organised as follows. In Section |2] 
we introduce the diagrammatic method which we use for 
the investigation of the single-band model. The class of 
diagrams which requires a special treatment due to their 
long-range contributions is discussed in Section [3] In Sec¬ 
tion 0] we show numerical results and focus, in particu¬ 
lar, on the convergence of our diagrammatic scheme. Our 
presentation is closed by a Summary and Outlook in Sec- 
tion|5] Some technical parts of the presentation are referred 
to three appendices 

2 Model and Method We investigate the single-band 
Hubbard model 


H = Ho + uY,di, (2) 

i 

in two dimensions, where 

Ho = J2 > di = . (3) 

ij.cr 

Here, i = (^ 1 ,^ 2 ) denotes one of the L sites on a square 
lattice, and a =t, i- The properties of this model will be 
studied in the thermodynamic limit T 00 by means of 
the variational wave functions 

Itf'G) = ^g|«'o) =nAG|tf'o) , (4) 

i 

first introduced by Gutzwiller l|4|, where |ifb) is a (nor¬ 
malised) single-particle product state and the local ‘Gutz¬ 
willer correlator’ is defined by 

Pi =^Ar|T)ii(r| . (5) 

r 

It contains the variational parameters A r for the four local 
states 

|P)ie{|0)i,lt)i,U)i,in)i} (6) 

for the empty, singly, or doubly occupied site i. Note that in 
Eq. (|5]l we have already assumed a translationally invariant 
ground state which allows us to work with parameters Xp 
that do not depend on the lattice site i. 

The single particle state jPg) is also a variational ob¬ 
ject and may be chosen as the ground state of an effective 
single-particle Hamiltonian, 


iJ,o- 


+ (7) 

i/j 


The effective hopping and pairing parameters f®? and Z\?? 
can then be considered as variational parameters which de¬ 
termine jPo)- Note that, in the main part of this work, we 


will consider superconducting ground states with d-wave 
symmetry for which the local pairing amplitude vanishes. 



0 — 



( 8 ) 


Here we introduced the notation (.. .)o,g for expectation 
values with respect to jPg) and |Pg)- The case of a finite 
local pairing ([8]) is discussed in AppendixlAl 

2.1 Diagrammatic expansion We need to evaluate 
the expectation value of the Hamiltonian (|2]), 


■^G = ^ fij 

hj,<y 


(’^G|c|',^Cj,^|pG) 


Y ('^gI'^'g) 


(9) 

with respect to our Gutzwiller wave function (IHi. As first 
shown in Ref. m, we can develop an efficient diagram¬ 
matic scheme for this evaluation if we demand that 


p/Pj = p2 = l + , (10) 


where 


?HF _ -HF-HF -nt _ - _ ^ 

Ui — ? ^l,cr — ^l,cr ^0 


HF _ , 


( 11 ) 


and ng = {n\^a)o = N/{2L). Equation (fTOl i determines 
three of the four parameters Ar as well as the coefficient x. 
In this way, we are left with only one variational parameter. 
Eor instance, we may express the parameters Ar by the 
coefficient x, 

Xl = l+x{l-nof , (12) 

A^ = 1 - xno{l - no) , (13) 

A 0 = 1 -f xuq . (14) 

For the calculation of Q we need to evaluate three 
power series in x, 

00 ^ 

k—0 

00 ^ 

md,m = a^E|[ (i6) 

k—O 
00 ^ 

{d'Glcl^C. JWo) = Efr (17) 

k—0 

where we used Eq. (fTOl) and introduced the notation 

dfZ.A.=d?^---df^ > dT^ = l, (18) 

^ AcgA . (19) 

The primes in Eqs. (fTSt - ifTTI i indicate the summation re¬ 
strictions 


Ipy^lp/, Ipy^iJ Vp,p'. (20) 

The expectation values (.. .jg in (fTSt - ifTTl i can be eval¬ 
uated by means of Wick’s theorem ifTOll . In the resulting 
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diagrammatic expansion, the fcth-order terms correspond 
to diagrams with k ‘internal’ vertices on sites li,... ,1^, 
one (two) ‘external’ vertices on site i (i and j) and lines 

= iiA'Jo , (21) 

= (cl,t4.i)o = (ci'4Ci,t)S (22) 

connecting these vertices. By construction, we eliminated 
all diagrams with local ‘Hartree bubbles’ at internal ver¬ 
tices, i.e., diagrams with lines that leave and enter the same 
internal vertex. To achieve the same for the external ver¬ 
tices in (fT6l l. (fT7l i we rewrite the corresponding operators 
as 


di = (1 - xdo)d^^ + no(n^f + n-Pf) + doP? , 

(23) 

S.IT "A,IT ^ ’ 

(24) 

where we introduced 


do = P-g , 

(25) 

q = Ai(Adng -1- A0(1 - no)) , 

(26) 

a = Ai(Ad — Ag) , 

(27) 


and f =1, i =t. When inserted into (fT^ . the last term 
in (l23t combines to A^(io('f'G|!^G) so that it does not have 
to be evaluated diagrammatically. 

As a result, we obtain diagrammatic sums with no 
Hartree bubbles at any vertex. This allows us to replace 
the lines (1211 1 by 

Aa' = - Siyno ■ (28) 

As demonstrated in Ref. Ga, the elimination of Hartree 
bubbles has significant consequences for the convergence 
and, hence, accuracy of our diagrammatic expansion. Due 
to our assumption of d-wave superconductivity, we do 
not have to eliminate ‘anomalous’ Hartree bubbles of the 
form ([^i. In AppendixlAlwe explain how our diagrammatic 
method can be generalised if Eq. (|8ll is not fulfilled. 

As the final analytical step of our derivation, we ap¬ 
ply the linked-cluster theorem lfT9l . The norm (fTSl l can¬ 
cels the disconnected diagrams in the two numerators (fTbl l 
and (fTTl i. Note that for the application of this theorem, we 
first need to lift the summation restrictions in Eqs. ([I5]l- 
dnii. This can be done, however, without generating addi¬ 
tional terms, as we explain in AppendixiBl 

Eor a translationally invariant system, the remaining 
task is to evaluate the diagrammatic sums 

OO 

^ = T.^S{k) (29) 

k=0 


with 


g ^ 2^_())i(3) ji(3)!(1) 2 ^^( 3 ) i ( 3 )|^ 


(30) 


and 

imm^k)^ Y: (31) 

ya)[(3)],(i)[(3)]|.^^ 

= Y ■ (32) 


Here, (•. • )o indicates that only connected diagrams are to 
be kept. Note that in the evaluation of these diagrams, i.e., 
after the application of the link-cluster theorem, one must 
not use any summation restrictions as in (fTbl l and (fTTl i. see 
AppendixiBl 

The structure of the variational ground-state energy 
functional is the same as in the paramagnetic case IH and 
given by 


{H)G = EG{\<I'o),x)=L{E'^^^ + Ud) (33) 


where 


ij 

+L[/A^((l - + do) . (34) 


This energy has to be minimised with respect to It^b) and x 
where |!f'o) enters the energy expression solely through the 
lines (I22I 1. (l28l l and through no- Note that in the presence 
of superconductivity, the particle number per lattice site 


nG = (ni,a)G = >i{do + - xdo) + 2no/^^^) (35) 

-l-A^ (m^ -I- (1 - 2no) - (1 + xrn()) , 


with nil = ^0(1 ~ iT'o) is not the same as ng. Physically, 
however, the value of riQ rather than no should be fixed 
in the minimisation of the energy. Therefore, we minimise 
the grand-canonical potential 


E = Eg — 2fiGnGL (36) 


with respect to Pfi>, S\,y, no, and x where the chemical 
potential pc allows us to vary the correlated particle num¬ 
ber. The minimisation with respect to and uq 

leads to the effective single-particle equation for 


Hf\Eo) = Eo\Eo) 


(37) 


with a Hamiltonian TTg® as introduced in (|7]i and parame¬ 
ters 




IJ 


dPi,i 


,eff ^ dE{\<Eo),x) 

dno 

ff dE{\Eo),x) 


(38) 

(39) 

(40) 
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Figure 1 Lowest order diagram of a) 7^^^, b) 7^^^ and ex¬ 
amples for long-range diagrams in c) and d) 


Equations (l34l i- (l40l i need to be solved self-consistently, 
together with the minimisation condition 

(41) 

Numerically, this can been achieved by the following iter¬ 
ative procedure; 

(i) Chose an initial value for ITTq) = \^o)- 

(ii) Determine the variational parameter a:inin which min¬ 
imises Eg for a fixed |t7b) = |!7o)- 

(iii) Determine the parameters (El-® and the corre¬ 
sponding Hamiltonian 77g® for x = ccmin- 

(iv) Determine the ground state j'T'g) of . 

(v) If ItT'g) ~ Jt^o) terminate the algorithm. Otherwise, set 
I'T'g) = I'T'g) and go back to point (ii). 

Note that, numerically, it is usually necessary to introduce 
some form of ‘damping’ in the calculation of the Hamilto¬ 
nian 77g® in step iii); If 77®® has been used in the previous 
iteration and 77^® is the Hamiltonian derived from (l38l l- 
(l40t then one continues the procedure in iv) with 

Hf = Hf + P{Hf - Hf) . (42) 

Working with a parameter /3 < 1 ensures the convergence 
of our algorithm. 

2.2 Calculation of diagrams To carry out the min¬ 
imisation, as described in the previous section, we need to 
calculate the diagrams (El and their derivatives with re¬ 
spect to lines up to a certain order in x. For example, the 
first-order diagram 7^^^ of 7^^^) is shown in Fig.[T^). Here, 


the lines can be normal or anomalous (S'],!!). With 

four normal lines, e.g., we have to evaluate 


7 ^) = 


pt pt P'l pi 


(43) 


= ^ ^ ?^k,t»^k',t»^k",4,?^k+k'+k",4, - nl (44) 

k.k'.k" 


where we have introduced the momentum-space distribu¬ 
tion rik.CT = (cjj. o-)o- Obviously, the real space eval¬ 
uation of the diagram is numerically much easier because 
only one summation (over li) has to be carried out. More¬ 
over, the lines Pij in real space vanish like — j| 

while the number of neighbours of this distance is ~ | i — j | ■ 
Hence, the real-space contributions of 7j^^ fall off rapidly 
and we can restrict the summation over Zi to a limited num¬ 
ber of nearest neighbours of i. This ‘locality’ of diagrams 
in real space is a key ingredient in our numerical imple¬ 
mentation since it allows us to calculate diagrams up to 
relatively large orders in x. 

Unfortunately, not all diagrams are as local as 7^^^. In 
particular, all diagrams in I^‘^\ contain ‘long-range contri¬ 
butions’, e.g., the joint sum over li, I2 in Fig. [T]?). In the 
paramagnetic case, there exists a relationship between 7^^^ 
and 7*^^^) which can be used to circumvent the long-range 
contributions in 7^^^ see Ref. na. No such relationship, 
however, can be used for superconducting states. More¬ 
over, even for paramagnetic states some diagrams have 
long-range contributions, e.g., the 7*^^^ and dia¬ 

grams shown in Fig.[T] We may identify the diagrams with 
long-range contributions by a topological analysis, as we 
shall explain in the following Section E] The real-space 
summations which belong to such long-range diagrams can 
then be evaluated analytically, see below. 


3 Long-range diagrams As explained in Sec¬ 
tion 12.21 some diagrams are not localised and require a 
special treatment in our real-space evaluation. Topologi¬ 
cally, there are two types of diagrams which we need to 
consider up to the 4-th order (of internal vertices) in x. 
They are displayed in Figs. |2] Diagrams of type I can be 
split into two disconnected diagrams Di and D 2 by cutting 
two lines, where the external vertices (i, j) belong to Di. 
Examples for such diagrams are shown in Figs. [TJ))-d). In 
a similar way we define diagrams of type II as those which 
can be split into three disconnected diagrams by cutting 
three (single) lines. An example for this type is the 
diagram shown in Fig. |3] It illustrates that type II diagram 
can only appear if there are at least four internal vertices. 
Note that more complicated long-range diagrams require 
the inclusion of more than five internal vertices. 

For the evaluation of the long-range diagrams in 
Figs. El one can carry out the sums over 1 (type I) or li, I2 
(type II) analytically. We will consider the paramagnetic 
and the superconducting case separately in the following 
two sections. 
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In a long-range diagram of type II we need to evaluate 


Figure 2 Long-range diagrams of a) type I, b) type II. 



Figure 3 A type II diagram in . 


3.1 The paramagnetic case For the evaluation of 
the diagram in Fig.|2^), we need to calculate 

Di='£ E KA+m'.m (45) 

m,m' 1 

where D^, are assumed to be localised diagrams, i.e., 
the sums over m, m' can be restricted to a shell around 0. 
The sum over 1 yields 

D{m, m') = E ^0,lA+m',m 

1 

= ^ ' (A.l ~ ~ 

1 

= (1 - 2no)P^ ^, + i5m,m'no • (46) 

Here we have used 

E^a^l+m',i+m=^^m,m' > (47) 

1 

which holds because, after Fourier transformation, in the 
paramagnetic case we can use ^ = nk.o- in momentum- 
space. A long range-diagram of type I is therefore given 
as 

77/= E (48) 

m.m' 

Note that the diagram may contain additional long- 
range elements as, e.g., in the diagram in Fig. H] In 
such a case, Eqs. (I46l l- (l48l l have to be applied consecu¬ 
tively. 


Dh 


where 


E ^Li,0^0,in2^m3,0 i:>(mi,m2,m3) 

mi.m 2 ,m 3 

(49) 


ll,l2 

(50) 

The sums over li, I 2 can again be calculated exactly. This 
leads to 


i9(mi,m2,m3) 


(1 3no-f 3nQ)P^2 

JT 3 

^7712 ,mi+m3 ^0 ■ 


(51) 


3.2 The superconducting case In the supercon¬ 
ducting case, the single (red) lines in Figs.|2]can be normal 
or anomalous. We therefore introduce the abbreviations 


Ej = Ki ’ 

(52) 

= *S'i,j , 

(53) 

-^Ot _ -\rOt C \rOL 

(54) 

and the corresponding Eourier transforms 


III 

(55) 

76k = (4,tC-k,4,)o ■ 

(56) 

Eor long-range diagrams of type I we then have to evaluate 


(57) 

m.m' 



(58) 


1 


where g {1,2} characterises the two single lines in 
Fig. 12^). Note that Xq q = 0 for our d-wave states. With a 
transformation to momentum space we find 

m') = - xi^XZ,,^, 

(59) 

with 

(60) 

k 


Note that, in the paramagnetic case, we have = 

m' ’ Tfo 0 = ^ 0 , such that Eq. (l46l l is recovered. 
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Figure 5 The 7-th order Taylor expansion of (l65l l as a func¬ 
tion of doping 6 for U/\t\ = 10 and different values of the 
cutoff parameter Tc = 1 — 10 (calculated without LRDE). 



Figure 6 Difference between the exact value of the dia¬ 
gram (see Fig.lTJ))) as a function of doping S for sev¬ 
eral values of Tc with and without LRDE (solid and dashed 
lines, respectively). 


The evaluation of type II diagrams leads to 




where 


mi ,m 2 ,m 3 

012,013 m 3 ) 


( 61 ) 


4 Results In our real-space evaluation of diagrams 
there are two main approximations that are needed to make 
the problem numerically treatable: 

i) The lines Pfy, S\^\' which enter the energy functional 
can have an arbitrary ‘length’ 


(mi, m 2 , m 3 ) 

I 1 .I 2 


The momentum-space evaluation for (l62l i yields 

-[KhYZ:Z+m3 + p^^^-] 

+ ["^o,o"^o,o"^m2,mi-i-m3 + perm.] 




VOLl 

i+ms^O.O^O.O^O.O 


(62) 


(63) 


where ‘-fperm.’ denotes the three cyclical permutations of 
(ai, a 2 , 03 ) in the respective functions and 


1 


^011,0:2,0:3 — J |^03 gi! 

m.m r k k k 


k(m—m^) 


(64) 


|i-i'| = ^(Zi-z ()2 + (/ 2 -;') 2 . 

To keep the number of real-space contributions finite 
we need to introduce some ‘cutoff’ r^, i.e., the assump¬ 
tion that P^y = = 0 for |1 — I'p > r^. The cut¬ 

off leads to numerical errors in particular for the long- 
range diagrams discussed in Section[3] However, as we 
will demonstrate in the following section ItTI these er¬ 
rors are negligible if we employ the long-range dia¬ 
gram evaluation (LRDE) technique introduced in Sec¬ 
tion [3 

ii) The number of diagrams grows exponentially with the 
number of internal vertices (index k in Eqs. OTl i. 
Therefore, the diagrammatic expansion must be termi¬ 
nated at some finite value of k. As we will discuss in 
section 10 a better truncation parameter is the total 
number of lines in a diagram. 


Note that in the superconducting case it is not possible 
to write the long-range contributions in terms of lines ( PA, 
j) as it is possible in the paramagnetic case, see Eqs. (l46l l 

and (l50l) . Instead, there appear the new objects Y^‘^, and 
whose calculation, however, is numerically be¬ 
nign because |m — m'j can be assumed to be small. Still, 
due to the appearance of these new objects, we need to 
reconsider the minimisation of our energy functional with 
respect to jt^o). This problem is discussed in ApnendixICl 


In all subsequent results, we have worked with a single¬ 
band Hamiltonian Hq that contains nearest and next- 
nearest neighbour hopping of f = —0.35eV and t'/t = 
—0.25, respectively. These values are generally assumed 
to describe the situation in the Cuprates. 

4.1 Line cutoff To analyse the role of a finite cutoff 
length Tc we first consider the expression for the correlated 
particle number per lattice site 

nG-no = [l + a;no(l - + x{l - 2no)P'^'> (65) 
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Figure 7 Condensation energy in Kelvin (upper panel) and 
correlated gap in the superconducting state (lower panel) as 
a function of doping 6 for U/\t\ = 10 and different values 
of Ic- 

which results from Eqs. (fTSl i. (fOT l. (l35l l. The correla¬ 
tion operator Pq in © commutes with the operator 
N = cr which counts the total number of elec¬ 
trons. Since, in the paramagnetic case, \Po) is an eigenstate 
of N, we have {N)g = {N)o and therefore no — n-o = 0 
in our translationally invariant system. 

If we consider the r.h.s. of dbSl) as a power series in x 
each coefficient of this expansion has to be exactly zero. 
Numerically, however, this is not the case because of our 
cutoff parameter Tc < oo. In Fig. |5] we plot the 7-th or¬ 
der Taylor expansion of dbST l as a function of doping for 
U/\t\ = 10 and different values of r^- 

As we can see from this figure, an increase of Tc im¬ 
proves the results only slowly. Unfortunately, for higher- 
order diagrams it would not be feasible to work with val¬ 
ues of Tc significantly larger than 10. Therefore, the LRDE 
is essential to eliminate the error that stems from the finite 
cutoff Xc- In fact, the LRDE ensures that the 7-th order ex¬ 
pansion in Fig.|5]is exactly zero for all cutoff parameters 
Tc- This perfect agreement, however, is only due to the fact 
that the remaining errors in the calculation of and 
cancel each other exactly in (1651) . To gauge the remaining 
real-space cutoff error, we display in Fig.|^the value of the 
diagram (see Fig. [TJi)) relative to its exact value as a 
function of doping for several values of Xc with and with¬ 
out LRDE. Note that these values are independent of x (cf. 
Eq. OTI) ) and only depend on the lines which we take from 
the ground state |!fh) of the bare single particle Hamilto¬ 
nian Hq. Again, we can see from Fig.|6]that bringing down 
the numerical error by increasing Xc is not working well 
without the use of the LRDE: even for Xc = 2 (i.e., only 
nearest and next-nearest neighbour lines) the results with 
LRDE are more accurate than those for Tc = 16 and with¬ 
out LRDE. 



Figure 8 Difference nc — n-o in the paramagnetic phase 
as a function of doping S for U/\t\ = 10 and several val¬ 
ues of Ic, with and without LRDE (solid and dashed lines, 
respectively). 

4.2 Line number truncation The natural expansion 
parameter appears to be the number k of internal vertices 
in Eqs. (l29]l- (l3Tl i. In fact, in our previous works on Pomer- 
anchuk phases M and superconductivity Ea we have in¬ 
vestigated the convergence of results as a function of the 
diagrammatic order k. However, the topological complex¬ 
ity of a diagram is more related to the number of lines in a 


diagram which is given as 


Nk = l + 2k 

for , 

(66) 

Nk = 2 + 2k 

for^ 

(67) 

Nk = 3 + 2k 

for T(3)(3)[fc] . 

(68) 


As we have found already in a study on the t-J model ifTSll 
it is more useful to include all diagrams up to a certain 
maximum number Ic of lines. This means that in results 
with Ic = 15 there are some diagrams included that have 
k = 7 internal vertices. As an example for the conver¬ 
gence with respect to Ic, we show in Fig. |7] the conden¬ 
sation energy (energy difference between superconducting 
and paramagnetic ground state) and the ‘correlated gap’ 
Aq = (ci_-|'Cj. 4 ,)G (for nearest neighbours i, j) in the su¬ 
perconducting state for U/\t\ = 10 and as a function of 
doping for different values of Ic- All these data have been 
calculated using the LRDE (without LRDE we obtain qual¬ 
itatively similar behavior, with the value of the conden¬ 
sation energy slightly increased, cf. also Fig. |9l). As ob¬ 
served in previous studies, convergence of the condensa¬ 
tion energy is reached for doping S > 0.1. For smaller 
doping values the convergence is less satisfactory as com¬ 
pared to that of most other observables, e.g., the correlated 
gap. Figure |7] shows that the results for the latter have con¬ 
verged already for Ic = 11. Since the condensation energy 
is largely increasing as a function of Ic, the stability of a su¬ 
perconducting state is very likely in the inaccessible limit 
Ic oo. Note that the ground state energy of the para- 
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magnetic phase (not plotted) is practically converged for 
Ic = 15 (the differences in this energy between the results 
for Ic = 15 and Ic = 11,13 are below IK). Therefore, the 
error of the condensation energy comes mostly from the 
superconducting phase ground state energy. In the follow¬ 
ing analysis we work with 4 = 15 and Tc = 10, unless 
stated otherwise. 

The importance of the LRDE is illustrated, once more, 
in Fig. [8] where we display uq — ng in the paramagnetic 
phase as a function of doping for U/\t\ = 10 and several 
values of Ic- The solid (dashed) lines show the results with 
(without) the LRDE. Obviously, the error that appears in 
the data without LRDE is so large, that there is hardly any 
improvement of the results if we increase 4- In contrast, if 
we use the LRDE, nc — no goes exponentially to zero if 
we increase Ic- Note that the remaining error stems from 
the fact, that contributions from the highest order in x are 
not exactly cancelled as they were in the Taylor-series ex¬ 
pansion discussed in the previous section I tTI 

In Fig. |9] we display the differences between energies 
that were calculated with and without the LRDE as a func¬ 
tion of doping and for U/\t\ = 10. The differences shown 
in this figure are those for the kinetic, the potential and the 
total energy in the superconducting phase, as well as the 
condensation energy. As can be seen from this graph, the 
changes of the energies due to the LRDE are not negligi¬ 
ble, in particular, for a doping of less than 0.1. The total en¬ 
ergy is lowered by up to 23K if we use LRDE. The result¬ 
ing change in the condensation energy, however, is smaller 
(less than lOK). The same holds for other observables in 
the superconducting state. Therefore, the main results on 
superconductivity which we have published in our previous 
work Ea remain unchanged by the new LRDE scheme. 
However, the relatively large changes of the kinetic and the 
potential energy indicate that for other systems or states, 
the LRDE may alter physical properties more visibly. 

In Fig. [To] we show the difference between the ki¬ 
netic energy of the paramagnetic and the superconduct¬ 
ing phases. This property is negative for a conventional 
(BCS-type) superconductor, as pairing induces ’smearing’ 
of the single-particle distribution around the Fermi surface, 
which increases the kinetic part of the total energy. We ob¬ 
serve such behaviour for U ;< 12|f| at all doping values. 
For larger values of the Coulomb interaction (U > 13|f|) 
the kinetic energy becomes lower in the superconducting 
phase. This is an unconventional behaviour coming from 
the fact that the bandstructure of the effective Hamiltonian 
changes upon condensation, which can dominate over the 
mentioned effect of an increased kinetic energy. The phe¬ 
nomenon of kinetic-energy driven superconductivity has 
also been observed experimentally for the cuprates II20II211 
[22ll2^ . However, the experimental trend is different in the 
sense that there is a transition close to optimal doping from 
BCS-type behaviour (for large doping values) to kinetic- 
energy driven superconductivity (for small doping values). 
In our previous calculations lITSll we obtained similar be- 



Figure 9 Differences between kinetic (Eikin), potential 
(Epot), total (Eg), and condensation (AE) energy calcu¬ 
lated with and without the LRDE as a function of doping S 
and for t//|f| = 10. 

haviour for U > 13|t|, however the kinetic energy increase 
was always very small. The present results are more ac¬ 
curate and we observe that the kinetic energy change for 
overdoped systems is positive for U > 13|f|. Close to 
the critical doping the kinetic energy change is very small 
(of the order of IK). It might be below the accuracy of 
our method to determine whether it is positive or negative 
in this regime (by changing Ic we could observe both for 
Ic = 11,13,15). The lack of transition between the two 
regimes in Gutzwiller Wave Function has been remedied 
in the variational Monte Carlo method llTlISl by including 
in the projection also an additional Jastrow factor. This Jas- 
trow factor is motivated by the form of a strong coupling 
expansion used to derive the t-J model from the Hub¬ 
bard model ll24l . By including similar terms in the present 
method it should be possible to observe a behaviour con¬ 
sistent with the experimental data. Work along this line is 
planned in the future. 

5 Summary and Outlook In summary, we have 
given a comprehensive derivation of a diagrammatic varia¬ 
tional method for the evaluation of superconducting states 
in two-dimensional Hubbard models. Since most diagrams 
in our scheme are rather localised in real space we are 
able to evaluate them up to relatively large orders in the 
expansion parameter x- For those diagrams which are not 
localised, we developed a resummation method that prac¬ 
tically eliminates the numerical error in their real space 
evalution (e.g., we estimated this error to be of the order of 
1 K for the condensation energy). The remaining error of 
our method comes from the cutoff in the expansion in x- 
We have analyzed convergence of the condensation energy 
and correlated gap as a function of this cutoff. We have 
also studied the kinetic energy change upon condensation, 
a property that can be related to the experiment. 
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Figure 10 Kinetic energy change upon condensation 
Z\i?kin as a function of doping 6. 


Our diagrammatic method is rather general and can 
be applied to various systems and situations. For the two- 
dimensional Hubbard model, it is still an open question 
whether a Pomeranchuk and a superconducting state are 
competing or coexisting near half filling. Also, the com¬ 
petition or coexistence of these two phases with antiferro¬ 
magnetic order Il25l has not yet been studied. 

It is further possible to apply the method to more com¬ 
plicated model systems, in particular to those, in which 
methods based on the Gutzwiller Approximation have 
provided valuable insights, e.g., periodic Anderson mod¬ 
els 026II27I . multi-layer Hubbard models, or multi-band 
Hubbard models 028II29I . Work in all these directions is in 
progress. Also the study of non-local interactions and/or 
correlations should be feasible in the future. 
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A Treatment of local pairing In this Appendix, we 
explain how our diagrammatic formalism can be applied if 
the local pairing condition (O is not fulfilled (e.g. for su¬ 
perconducting phase with s-wave symmetry component). 

The basic idea remains the same as before, i.e, we aim 
to write P/Pj in the form (0 with an operator that 
ensures the vanishing of all Hartree bubbles at internal ver¬ 
tices. Now, however, we need to cancel normal as well as 
anomalous Hartree bubbles. This requires the use of the 


more general local correlation operator 


Pi = ^Ar|P)ii(P|+AB(|d)ii(0|+h.c.) . (69) 

r 

Here we have already assumed that 

Aq = (4i)o (with Ai = ci,|Ci,t) (70) 

is real which allows us to also work with a real parame¬ 
ter As. The following considerations can be readily gener¬ 
alised in the case of a complex amplitude Z\o- 

It will be useful to generalise the Hartree Eock opera¬ 
tors (fTTI) . Let 

Oi = ai.i ■ • ■ ai.n (71) 

be an arbitrary operator on site 1 where dpi may be creation 
or annihilation operators. Then, we want 

OjHF = o, _ [Oj]HF (72) 


to create the same diagrams as Oi apart from those with 
Hartree bubbles at site 1. This is achieved if we define 
recursively as 

[dpi... dp„]^^ = (dpi... dp„)Q (73) 


{7l.- -.7n}=0 y \f=l / 

xlfi-ir) 

\e=i / 0 

with 

/s({ 7 *}) = ^ 7 ^ • 


n 

HF 'j 


1 

7=1 

i 


(74) 


The prime in (l73t indicates that 


2 < y] 7^ < n - 2 (75) 


has to be even (odd) if n is even (odd). Due to this summa¬ 
tion restriction we find 

[dpidpj]^^ = (dpidpj)o , (76) 

[dp.]“^ = 0 . (77) 

Note that this recursive definition is quite general and cov¬ 
ers systems with an arbitrary number of orbitals and local 
density matrices (di 7 dij)o. In the case of our single-band 
model with local pairing (iTOl i. we find, e.g., 

jpP = d\ - [cppcppcp^cpp]^^ (78) 

= d\- no(npp -f npp) - Ao{Ai + A}) + nl + Al 

for the Hartree Eock operator in Eq. (ITOl i. Another exam¬ 
ple which will be relevant for the evaluation of hopping 
expectation values is 

[cj ^cppcpp]^^ = nocpp -f /AoCpt ' 
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Together with (fTsT l. the operator equation (fTOl i. leads to 


= 1 + (1 “ no^x + A^x , (80) 

Xl = I - nQ{l - no)x + Alx , (81) 

A0 + Ag = 1 + ttqX + Ai^x , (82) 

Xs^Xd + Ag) = —xZ\o . (83) 


This set of equations determines, like in the case without 
local pairing, all parameters Xr, Xb as a function of x. 

To calculate the expectation value of a local doble oc¬ 
cupancy, we use the relation 

P^diP, = {Xl + Xl)<^^ (84) 

+ [(A2+A|)no-A|](nPF + nPf) 

T[('^d T Xb)Ao — ABAd](^i^^ -f h.c.) + do 

with 

- Z\o (85) 

and 

do = (A^ + X^){no + 4\q) + Ag(l — 2no) -f 2XBXdAQ . 

( 86 ) 

This leads to 

{di)G = do + 2[iXl + Xl)no - (87) 

+('^d + ~ xdo)!^'^^ + 2[(A^ -f }?b)Ao — AsAdJ/p^ 

where we have introduced the (anomalous) diagrammatic 
sum 

oo ^ 

E (88) 

k=0 li,...,lfc 

For the evaluation of a hopping expectation value we 
expand (fT^ , as in Eq. (l24l) . 


+a^Ci gni^a - (89) 

Cucr = “ [4, 

+a(cl-ni^cy - ( 90 ) 

with 

q = Ai[Adno -f A 0(1 — no) + 2 XbAo] , (91) 

q = Ai[(Ad — A 0 )Z\o -f Ab(1 — 2no) + 2AbZ\o] , (92) 

a = Ai(Ad — A 0 ) , (93) 

a = -2AiAb . (94) 


Using the above equations, it is now a straightforward 
task to write down the somewhat lengthy expressions for 
{AA,<t)g and (c)^ct^)G. 


B The linked-cluster theorem In order to apply the 
linked-cluster theorem in Q we need to lift the summa¬ 
tion restrictions (l20l i in Eqs. (fT5l l- (fT7l i. As we shall explain 
in this Appendix, the summation restriction can be lifted 
without generating additional terms. 

We consider a diagram D in which the two operators 
Cl cr and Cl'cr from the two internal vertices 1 , T have a con¬ 
traction with two other operators di, 0 . 2 , respectively. The 
operators di need not to be specified, i.e., they could be 
creation or annihilation operators and may belong to inter¬ 
nal as well as external vertices. The diagram D results as 
one term in the evaluation of 

(cu .crCl.crdi, Q;2f9rest )o (95) 

by means of Wicks theorem where Orest contains all other 
operators which appear in D. Hence, we can write D as 

D = {c\,adi)o{c\> ^ad2)oD rest • (96) 

Another contribution from (l95T l, however, is 

P — (ci,rrCt2)o (ci'.rr tti)o79rest • (97) 

If 1 = 1', we have D + D' = 0, i.e., both diagrams can¬ 
cel each other. This explains why the summation restriction 
for the internal vertices can be lifted. The same arguments 
work if one of the two considered sites belongs to an ex¬ 
ternal vertex, since all four operators c^'l generate contrac¬ 
tions for each internal vertex. 

Note that after the application of the linked-cluster the¬ 
orem, one must not reintroduce the summation restrictions. 
Eor instance, the diagram D could be fully connected while 
in D' = the factor is cancelled by the norm. 

In such a case, the sum D -f P)'^^ which results after the 
application of the linked-cluster theorem will, in general, 
not be zero even for 1 = 1 '. 

C Minimisation with respect to 
C.1 Energy functionai without iong-range dia¬ 
grams If we ignore the non-locality of long-range dia¬ 
grams, our Lagrange functional depends on x and jif'o) 
where I'^o) enters the functional only via the elements of 
the single-particle density matrix 



Here we introduced 

P{ia),{jcr) = (Cjo-^j'( 8 )o (99) 

'^(ifT),(jCT) = ((^jCTCj(T)o ■ (100) 

Since p® is derived from a single-particle product state it 
has to obey the constraint p® • p® = p® . Like in the param¬ 
agnetic case we implement this constraint in the minimisa¬ 
tion with respect to p® by means of Lagrange parameters, 
see Refs. II31II32II . This leads to the self-consistent single¬ 
particle problem introduced in Eqs. (|7l), (|J7l) - (l40l) 
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C.2 Energy functional with long-range diagrams 

When we evaluate the long-range diagrams as described 
in Section |2 we obtain an energy functional that does not 
only depend on the elements of p® in real space, i.e., the 
lines X“j. It also depends on the ‘higher-order’ lines 
and Z?} which, on the other hand, are determined by 
the elements of p® in momentum space, i.e., the distribu¬ 
tions XQ. Therefore, the great canonical potential has the 
form 


^({Pi,jh{Sij};{xn) (101) 

Note that the lines Pi j, j enter the functional also via 
the evaluation of long-range diagrams, see Eqs. (l52li, (153, 
(l59l l. (l63l l. The minimisation of P with respect to Pij, 5'ij, 
and X^ leads to 

= Pois'd) (102) 


where 






i^j 


+ Z «k^k,<T + ( Z ^kCk,t^-k.i + h-C- ] (104) 


and 



‘4 5Pi.j ’ 

Z3- = 

‘J ’ 

(105) 

dJF 



~ dxl ’ 

A 

dXl' 

(106) 

Equation dlOlli then yields 



£k = Z(^k“+<’')^k (107) 

ct 

I \ ^ / l,a,a' I a,1,0: , o,o , 1 \ -v^a-v^o' 

+ Z^(^k +^^1. +^k ) 

o,o' 

^k = Z(^k“+<’')^k (108) 

O 

I \ ^ / 2,0,o' I 0,2,o' I o,o', 2 \-v^o' 

+ Z^(^k +^k +^k JT^kT^k ■ 

o,o' 

where 

1 ^ . ^ 

(109) 


^ c‘‘‘0~o ^ P 

C T r ~^-^ ^ r \ r \' ’ 


0 ‘\ 7 ' 0 ,« 


-I ___ f ) 

1,02,03 _ ^ ik(.i-i) ^ j: 

— T / 0701,02,03^ • 


^ 2 '“!’*^ 2,03 * 


( 110 ) 
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